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ABSTRACT 

This paper presents stability calculations made for 
a shock-free supersonic jet using the model based on 
parabolized stability equations. In this analysis the 
large-scale structures, which play a dominant role in 
the mixing as well as the noise radiated, are modeled 
as instability waves. This model takes into considera- 
tion non-parallel flow effects and also nonlinear interac- 
tion of the instability waves. The stability calculations 
have been performed for different frequencies and mode 
numbers over a range of jet operating temperatures. 
Comparisons are made, where appropriate, with the so- 
lutions to Rayleigh’s equation (linear, inviscid analysis 
with the assumption of parallel flow). The comparison 
of the solutions obtained using the two approaches show 
very good agreement. 

1. INTRODUCTION 

The aerospace engineering community has placed a 
major emphasis on the development of a high speed civil 
transport aircraft. The success of this effort hinges on 
several factors. One of the major problems that needs 
to be addressed is the environmental issue of the noise 
generated by aircraft at take-off, cutback, and approach 
engine power settings. Jet noise represents a significant 
percentage of the overall noise generated. Hence, the 
reduction of jet noise is of utmost importance to fu- 
ture development of a supersonic civil aircraft. This 
has led to a revival of interest in developing techniques 
to predict and ultimately suppress noise radiated by su- 
personic jets. 

Several approaches have been developed for the 
prediction of jet noise. One such approach is the direct 
numerical simulation of the full Navier-Stokes equations 
for the estimation of the radiated noise. This technique, 
although exact, is not practical because of the require- 
ments of the computational resources. This has led to 
development of models based on simplified governing 
equations. Mankbadi et al. 1 developed such a technique 
by using large-eddy simulation to model the flow field 


f Research Associate, Member AIAA 

* Senior Research Scientist, Associate Fellow AIAA 

* Eminent Professor, Associate Fellow AIAA 


and calculated the far-field noise by applying LighthilFs 
acoustic analogy. These techniques also require exten- 
sive computational power. Recently, there has been 
lot of focus on extending computational fluid dynam- 
ics techniques for aeroacoustic applications. The de- 
velopment of these techniques are not simple as the 
requirements on numerical dissipation, dispersion, etc. 
for acoustic analysis are very rigid, see Tam 2 . 

A different approach for the prediction of jet noise 
is based on combining theoretical and analytical meth- 
ods. It is well established that turbulent flows possess 
coherent large-scale structures. These structures con- 
trol the turbulent mixing process and also play a dom- 
inant role in the noise generation process. It has been 
successfully shown, in free shear layers such as plane 
mixing layers and jets, that the large-scale motions can 
be modeled as a random superposition of instability 
waves. The noise generated by these waves compare 
well with the experimental data. A good review of this 
approach is presented in Tam 3 . This model can be used 
to determine the two important characteristics of an in- 
stability wave, phase velocity and growth rate, which 
determine the mixing as well as the noise radiated. 

Most of the earlier theoretical studies of supersonic 
jets have considered linear, inviscid analysis and have 
assumed the mean flow to be parallel. Tam and Burton 4 
took into consideration the non-parallel flow effects by 
using the method of multiple-scales. The real physi- 
cal process of jets involve nonlinear interaction among 
all turbulence scales. The nonlinearities considered in 
this paper may be divided into nonlinear interaction be- 
tween the instability waves and the mean flow, and the 
nonlinear self interaction of the instability waves. The 
former leads to a change in the mean-flow distribution, 
which in turn will modify the axial development of the 
instability waves. This kind of nonlinear effect can be 
taken into account by using measured mean velocity 
data. However, this data is not always readily available 
for different flow conditions. The effects of nonlinear 
self-interaction of the instability waves have been ne- 
glected in all the earlier studies by assuming that the 
amplitude of the disturbances, which initiate the in- 
stability waves near the nozzle exit, are small. This 
assumption may not be valid in all cases. 

In this study, nonlinear stability analysis of super- 
sonic jets is considered. This analysis also accounts 
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for non-parallel flows. The nonlinear interaction of the 
waves are predicted by the solving the parabolized sta- 
bility equations (PSE). Bertoloti et al. 5 and Chang et 
al. 6 have applied PSE based technique for the stabil- 
ity analysis of incompressible and compressible bound- 
ary layers. The mean flow varies slowly in the axial 
direction enabling one to decompose the perturbation 
quantities into a rapidly varying wave-like part and a 
slowly varying shape function. The PSE can then be 
derived from the Navier-Stokes equations by applying 
parabolizing approximation to the shape function. The 
simplified governing equations can be solved by march- 
ing along the streamwise direction. This method is not 
computationally expensive. 

In this paper stability calculations of supersonic 
jets using the PSE are carried out. The calculations 
have been performed for different frequencies and mode 
numbers over a range of jet temperatures, from cold to 
hot. In section 2, the development of the PSE model 
is presented and in section 3 computational procedure 
is described. The numerical results obtained for differ- 
ent cases are presented in section 4. Finally, in section 
5, a summary of the results is given and also ways of 
improving the model and plans for future work are dis- 
cussed. 
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The governing equations (2.3) have to be parabo- 
lized to enable a marching solution in the streamwise 
direction. This is accomplished by decomposing the 
disturbances into a rapidly varying wave-like part and 
a slowly varying shape function. We assume that the 
disturbance vector <)> for an instability wave with a fre 
quency mu> and azimuthal mode number n can be ex 
pressed as 


2. ANALYSIS 

Consider the development of instability waves gen- 
erated in a shock-free jet. It is assumed that the insta- 
bility waves are governed by the compressible equations 
of motion, energy and state for a perfect gas. The flow 
variables are decomposed into a mean value and a per- 
turbation quantity 

f(r,9,x,t) = f(r,0,x) + f(r, $,x,t) (2-1) 


where / is any flow variable and (r, 0 , x) is the polar co- 
ordinate system with x in the axial direction. The gov- 
erning equations for the disturbances are obtained by 
substituting the form given by equation (2.1) in the gov- 
erning equations and subtracting the equations for the 
mean flow. The resulting equations can be expressed as 
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where <j) is the disturbance vector and is given by 0 = 

( p ti', t/, w\p f ) T . The coefficient matrices, I\ A, B y ...» Vq x 
and V$$ are composed of a linear part and a nonlinear 
part. The linear part (denoted by superscripts /) con- 
tains only mean flow quantities and the nonlinear part 
(denoted by superscripts n) contains disturbance quan- 
tities , i.e. T = + T n , etc. Equation (2.2) can then 
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where the wave-like part \mn is given by 
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Q mn is the axial wavenumber and is the shape 

function vector given by = (/>, u, 0, w,p) T . As 

the shape function is assumed to vary slowly with the 
streamwise direction x, the second derivative d 2 ^ /dx' 
is very small and is neglected. The governing equations 
for the shape function of a single mode (m,n), using 
the form given by equation (2.5), can then be written 
as 
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where 
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The nonlinear forcing function F mn can be evaluated 
from the Fourier series expansion of F n 


m=oo n=:oo 

F n (r,9,x,t)= *22 exp[i(n^-mwi)] 
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In equation (2.7), the elliptic effect associated with the 
wave part is absorbed in matrices D y A and C and does 
not contribute to the upstream influence. However, the 
term dp/dx allows upstream influence in subsonic re- 
gions of the flow. In order to make equation (2.7) truly 
parabolic, dp/dx is multipled by a constant given by 

£7 = ( i+(7-i)a/2. Mx < 1 (2.9) 

l 1 M x > 1 

where M T is the local Mach number. This parabolizing 
procedure yield solutions for the parabolized Navier- 
Stokes equations which compare well with those ob- 
tained by the full Navier-Stokes equations, see Rubin 


3. NUMERICAL SCHEME 

The parabolized stability equation (2.7) is solved 
numerically by finite-differencing the derivatives with 
respect to x and r. In the formulation given below, 
j represents the shape function with an axial location 
index i and radial location denoted by index j. The 
streamwise derivative at (i y j) is replaced by a second- 
order backward difference given by 

d^ 

— = (3*,,; - 4^-ij 4- »i-aj)/2A* (3.1) 

for all axial locations except for the starting plane where 
a first-order backward differencing is used. The radial 
derivatives are discretized using a fourth-order central 
differencing and are given by 

!fr “ 12Ar 1 
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The above finite differencing scheme is replaced by a 
second-order scheme for the grid point next to the bound- 
ary. The second-order scheme can be written as 
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The finite difference formulations given above are sub- 
stituted into equation (2.7). This coupled with the 
boundary conditions results in a block penta-diagonal 
system of equations (block size of 5 x 5) at each axial 
location for the unknown shape functions, 

The boundary conditions have to be prescribed for 
the eigenfunctions at the jet centerline and at a radial 
location of infinity. As the governing equations become 
singular at r = 0, L’HospitaPs rule is applied at the 


centerline to obtain a new set of governing equations. 
The boundary conditions along the centerline takes a 
different form based on the azimuthal mode number. 
For example, for axisymmetric mode, i.e. n = 0, the 
boundary conditions at r =0 are given by 

u — 0 = 0, dw/dr = dp/dr — 0 (3.6) 


Here, the velocity components u y v y w are in the r, 6,x 
directions, respectively. The boundary condition on 
density is obtained from (3.6) and the discretized con- 
tinuity equation. At the outer boundary, far from the 
centerline, all the shape functions are set equal to zero. 
For the nonlinear analysis, the far-field boundary con- 
ditions can be directly applied for all modes except for 
modes with zero frequency (m = 0) and axisymmetric 
mode (n = 0). The case of m = 0 and n = 0 is referred 
to as the mean flow correction mode. For this mode, 
the condition u = 0 is replaced by du/dr = 0. 

The solution procedure is as follows: start with a 
known solution for <x mn and ^ m n at x — xq, assume 
a value for a mn at the next axial location, march the 
PSE to the next station by solving for ^ mn . The new 
value of a mn is determined from 


&mn, new 


C*mn, old 


i ran 

^ mn dx 


(3.7) 


and this process is repeated until the difference between 
the assumed and computed values of a mn is less than a 
given tolerance. The updating procedure given above is 
equivalent to normalizing the shape function such that 
dW/dx is zero at a particular radial location. As the 
shape function vector is a function of r and contains 
five dependent variables, there is no unique wavenumber 
as in the liner theory for parallel flows. However, as 
shown by Chang et al. 6 , the total growth rate appears 
to be very weakly dependent upon the normalization. 
Further details of the normalization and its effect on the 
computed wavenumber are given in Chang et al. 6 . In all 
the results presented here, we have chosen to compute 
the wavenumber based on streamwise velocity using 


new 


— a mn, old 


. f n °° u'dit/dx dr 


-where * denotes complex conjugate. 


(3.8) 


4. NUMERICAL RESULTS 

The preliminary results from this study are pre- 
sented here. The PSE has been used to analyze the 
linear stability of a supersonic jet issuing from a circu- 
lar nozzle operating at design Mach number of 2.0. It 
has been assumed that the axial evolution of the insta- 
bility weaves can be treated as an inviscid phenomena 
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and Fmn in equation (2.7) is zero. The results obtained 
for parallel flows using the PSE are compared with the 
solutions to equation governing linear, inviscid insta- 
bility wave, i.e., the compressible Rayleigh’s equation. 
The solutions to the Rayleigh’s equation (RES) are ob- 
tained using the method of Seiner, Bhat and Ponton 8 . 
For both the approaches, i.e., PSE and RES, the mean 
flow field is assumed to be known. It has been shown 
that the mean velocity profile can be closely approxi- 
mated by a half-Gaussian function given by 


W(x) 


W c (x) 

W c (x)exp[— ln2(r — h)/b 2 ] 


r < h 
r > h 


(4.1) 


where /*(x) is the radius of the potential core, b(x) is the 
half-width of the shear layer and W c (x) is the jet cen- 
terline velocity. For the present calculations, the values 
of these parameters are obtained from computational 
data of Viswanathan et al. 9 . 

Calculations have been performed for a range of 
frequencies and the modes considered are n = 0, the 
axisymmetric mode, n = 1 and 2, the helical modes. 
These calculations were carried out over a range of jet 
operating temperatures. However, for conciseness, the 
results for n = 1 mode are shown here for a cold jet. 
The two most important characteristics of an instability 
wave are the phase velocity and growth rate. The phase 
velocity is defined as as 2vffot r where / is the frequency 
and a r is the real part of the complex wavenumber a. 
The growth rate is defined as —a, where a* is the imag- 
inary part of the wavenumber. Figures 1 and 2 compare 
the phase velocity and growth rates obtained using the 
two approaches at two different axial locations. It can 
be seen that over the entire range of Strouhal number 
(St — 2 frj/Wj) considered, the two solutions are in 
very good agreement. Similar agreements are also ob- 
tained for n = 0 and 2. 



Strouhal Number 

Figure 1. Comparison of phase velocities 



Strouhal Number 

Figure 2. Comparison of growth rates 

Figures 3 and 4 show the comparison as a function cf 
downstream distance for axisymmetric mode (n = 0) 
at St = 0.54 and the helical mode (n = 1) at St == 
0.49. Once again, the agreement is excellent for both 
the phase velocity and growth rate. As mentioned be- 
fore, calculations have also been performed for hot jets 
obtaining good agreement between the solutions using 
the two approaches. 



Axial Distance, X/R, 

Figure 3. Axial variation of phase velocities 

The results presented so far have been obtained fcr 
parallel flows. We now consider the solution of linear 
PSE with non-parallel flow effects. As in any mark- 
ing scheme, an initial solution is needed to start the 
computations. To provide a good starting solution, the 
multiple-scales analysis can be used which accounts fcr 
non-parallel flow effects. It has been shown by Chang 
et al. 6 that using any incorrect initial conditions result 
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Figure 4. Axial variation of growth rates 

in a transient phenomena close to the jet exit. Further- 
more, it was observed that the PSE approach is capable 
of recovering the correct solution in the marching pro- 
cess provided the input initial conditions are close to 
the correct solution. In this study, the linear parallel 
stability equations are solved and the eigenvalue and 
eigenfunctions are used as the starting solutions. 

The linear PSE with non-parallel flow effects have 
been solved for several cases. Figures 5 and 6 show the 
results for two of these cases. Once again, calculations 
are carried out for a cold supersonic jet operating at its 
design Mach number of 2.0. In figure 5, the variation of 
the wavenumber (real part) and the growth rate with 
axial distance is presented for axisymmetric mode (n = 
0) at St = 0.54. The results from parallel stability the- 
ory are included for comparison. As can be seen, the 
wavenumber and the growth rates start off being less 
than those based on parallel theory. After some initial 
transients, both the wavenumber and growth rate are 
larger than those from parallel theory. This implies that 
the wave based on non-parallel flow has a lower phase 
velocity than the wave based on parallel flow. In ad- 
dition, the wave corresponding to the non-parallel flow 
gets damped further downstream than that obtained 
using parallel flow assumption. The results for helical 
mode (n = 1) at St = 0.49 are presented in figure 6. 
Similar trends as in axisymmetric mode are observed. 
The only difference is that the wavenumber based on 
parallel and non-parallel mean flow are more or less the 
same. 

The eigenfunction distribution for the cases dis- 
cussed above are presented next. Figure 7 shows the 
radial distributions of the eigenfunctions for the axisym- 
metric mode, St = 0.54 at X/Rj = 0. The eigenfunc- 
tions shown are the radial velocity component, u, axial 
velocity component, w, and pressure, p. For this case, 
the azimuthal velocity component, v , is zero. As can 




Axial Distance, X/R, 

Figure 6. Comparison of wavenumbers, n = 1, St = 0.49 

be seen, the axial velocity component is very large and 
peaks in the shear layer. These eigenfunctions further 
downstream, at X/Rj = 10, are presented in figure 
8. The amplitude of the pressure eigenfunction has in- 
creased slightly while the amplitude of the axial velocity 
component has decreased substantially. The axial loca- 
tion for peak noise emission, based on the experimental 
observations of Troutt and McLaughlin 10 and Seiner, 
et al. 11 , occurs near each wave’s neutral point. This 
suggests that the amplitude of pressure would increase 
with downstream distance until the wave has become 
damped. 

The eigenfunctions for n = 1, St = 0.49 at two 
axial locations, X/Rj = 0 and X/Rj = 10 are shown 
in figures 9 and 10, respectively. In this case, the az- 
imuthal velocity component is nonzero. It can be seen 
that the trends are similar to that of the axisymmetric 
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Eigenfunctions 



Radial Distance, r/R, Radial Distance, rf Rj 

Figure 7. Eigenfunctions for n=0, £^=0.54 at X/Rj=Q Figure 9. Eigenfunctions for n=l, £<=0.49 at X/Ri— 0 



Figure 8. Eigenfunctions for n— 0, £<=0.54 at X/Rj = 10 

mode. The only difference is that the amplitude of the 
azimuthal velocity component, t?, has increased. Con- 
verged solutions to the nonlinear PSE could not be ob- 
tained. Currently, this problem is being investigated. 
Some of the possible reasons for not getting converged 
solutions are discussed in the next section. 

5. DISCUSSIONS AND CONCLUSIONS 

In this paper, a model for linear and nonlinear sta- 
bility analysis of supersonic jets is developed using the 
parabolized stability equations. The linear PSE has 
been used to analyze the instability waves of a super- 
sonic jet operating at design conditions, i.e., shock-free. 
For parallel flows, it has been shown that the PSE ap- 
proach and the solutions to Rayleigh’s equation are in 
very good agreement over a range of frequencies and 
azimuthal mode numbers. The axial development of 


Figure 10. Eigenfunctions for n=l, £< = 0.49 at X/Rj= 10 

the instability waves also compared well. These calcu- 
lations also considered a range of jet operating temper- 
atures from cold to hot. The important characteristics 
of the instability waves, i.e., phase velocity and growth 
rate, are predicted well using the PSE. 

The linear PSE model has also been used to study 
the non-parallel flow effects. The computed solutions 
show that the growth rate is higher and the wave’s phase 
velocity is lower than that corresponding to the parallel 
flow. This expected trend is seen for instability waves 
of different frequencies and mode numbers. We have 
assumed that the dynamics of the axial development of 
the instability waves can be treated as an in viscid phe- 
nomena. At distances far downstream from the jet exit 
the instability wave would no longer grow in amplitude 
and is damped. At these locations the governing equa- 
tions become singular and the integration contour must 
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be deformed around the critical point, see Tam 7 * * * * 12 and 
Tam and Morris 13 . In this study, the calculations are 
terminated close to the neutral point of the wave. 

The results presented here are preliminary and is 
by no means complete. In order to improve this model, 
there are several issues that need to be investigated. 
First of all, the solutions obtained from linear parallel 
stability equations have been used as the initial con- 
ditions for the marching scheme. The sensitivity of 
the initial conditions to the converged solutions will be 
studied. Solutions obtained using the multiple-scales 
analysis, which accounts for non-parallel flow effects, 
can be used as the initial conditions for the eigenvalue 
and the eigenfunctions. Secondly, to circumvent the 
problem associated with the damped inviscid wave and 
to study the effects of Reynolds number, viscous terms 
will be included in the analysis. 

As mentioned before, converged solutions to the 
nonlinear PSE could not be obtained. It appears that 
the initial amplitudes of the various waves considered 
play a critical role in getting converged solutions. Fur- 
thermore, the selection of the waves, with nonzero ini- 
tial amplitudes, also influences the iterative process. 
Currently, the effects of the initial amplitudes and the 
choice of the waves on getting converged solutions are 
being investigated. The solutions obtained from non- 
linear PSE analysis, after all the issues related to con- 
vergence have been identified and resolved, will be vali- 
dated by comparing it to solutions computed from non- 
linear numerical simulation, for example, the results of 
Viswanathan and Sankar 14 . 
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